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We reformulate the theory of Schwarzschild black hole perturbations in terms of the metric per- 
turbation in the Lorenz gauge. In this formulation, each tensor-harmonic mode of the perturbation 
is constructed algebraically from 10 scalar functions, satisfying a set of 10 wavelike equations, which 
are decoupled at their principal parts. We solve these equations using numerical evolution in the 
time domain, for the case of a pointlike test particle set in a circular geodesic orbit around the black 
hole. Our code uses characteristic coordinates, and incorporates a constraint damping scheme. The 
axially-symmetric, odd-parity modes of the perturbation are obtained analytically. The approach 
developed here is especially advantageous in applications requiring knowledge of the local metric 
perturbation near a point particle; in particular, it offers a useful framework for calculations of the 
gravitational self force. 
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I. INTRODUCTION 

There has been much progress, over the past few years, on the problem of finding the "self force" (SF) correction 
to the geodesic motion of point particles in curved geometries [1,2]. This had strong motivation from the need to 
accurately model the orbital evolution in astrophysical binaries with extreme mass ratios — a key type of sources for 
the planned gravitational- wave detector LISA (the Laser Interferometer Space Antenna) [3]. Yet, despite the fact 
that a formal framework for calculations of the gravitational SF exists since quite a while ago [4-7], most actual 
computations so far have been restricted to toy models based on scalar or electromagnetic fields, with very little 
progress on the gravitational problem of relevance to LISA. What held progress up, mainly, is a certain gap within 
the standard glossary of techniques comprising black hole perturbation theory. Here we intend to supplement the 
necessary piece of perturbative technology needed to facilitate calculations of the gravitational SF. Although our 
motivation bears primarily on the SF problem, we expect the perturbation approach to be developed here to be useful 
in a wider class of problems in mathematical and numerical relativity. 

What has been hindering calculations of the gravitational SF can be described, in most simple terms, as follows. 
The term "gravitational SF" refers to the effective local force exerted on a point mass particle through interaction with 
its own gravitational field. (The notion of a "point particle", in our general-relativistic context, builds on the clear 
separation of length scales in the extreme-mass-ratio inspiral problem — see, e.g., the discussion in [1].) Technically, 
the SF appears as a correction term, quadratic in the particle's mass, in the geodesic equation of motion. To obtain 
the SF one first needs to address the subtle issue of "regularization" : telling the part of the particle's field whose 
back-reaction affects the motion, from a remaining, singular piece that can be safely "removed" without affecting the 
motion. A formal procedure for identifying the correct "regularizablc" piece (for mass particles in any vacuum curved 
spacetime) has been established in works by Mino et al. [5], Detweiler and whiting [6], and others. 

The above regularization procedure involves the introduction of a suitable local frame for the particle, and also 
the imposition of an appropriate gauge condition for the particle's field, the latter being treated as a perturbation 
over the fixed background of the massive black hole. Recall here that the metric perturbation (MP) is subject to a 
gauge freedom, associated with infinitesimal coordinate transformations. The full MP induced by the particle is gauge 
dependent, and so is the form of its singular, "regularizable" piece. This local piece is most conveniently constructed 
within the so-called "Lorenz gauge": an analogue of the familiar Lorenz gauge of electromagnetism, wherein the 
covariant divergence of the (trace-reversed) MP is taken to vanish [see Eq. (3) below]. The Lorenz gauge condition 
conforms with the isotropic form of the singularity very close to the particle (as viewed from a suitable local frame), 
and the Lorenz-gauge MP correctly reflects the "inverse-distance" behavior of this singularity. This local isotropy of 
the Lorenz gauge is an essential feature that should not be taken for granted: Expressed in a poorly chosen gauge, the 
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physical point singularity may take a complicated form, that could render analysis intractable. Indeed, it has been 
shown that with the standard gauge choices commonly made in studies of black hole perturbations (see below) the 
physical point singularity can take a distorted form, and may even cease to be isolated [8]. 

Once the "correct" singular piece of the MP is identified, the construction of the SF proceeds by removing this piece 
from the full (retarded) field of the particle, and then calculating the force exerted on the particle by the remaining 
finite field. (The full field is obtained, in principle, by solving the perturbation equations with an energy-momentum 
source term representing the point particle.) A practical way to subtract the singular piece from the full field is 
offered by the mode sum scheme [7], in which one first decomposes the full MP into multipolc harmonics (defined 
with respect to the background black hole geometry) and then carries out the subtraction mode by mode — thereby 
avoiding the need to deal with singular fields. The mode- by- mode subtraction scheme, like the original construction 
in [4,5], is formulated in the Lorenz gauge, and relies on having at hand the full Lorenz-gauge MP near the particle. 

Unfortunately, the standard toolkit of black hole perturbation theory does not include practical techniques or 
working tools for calculating the MP in the Lorenz gauge. Standard formulations of black hole perturbations employ 
other gauges, which are favored for their algebraic simplicity. A most commonly favored gauge choice for analyzing 
perturbations of spherically symmetric black hole spacetimes is the one introduced long ago by Regge and Wheeler 
(RW) [9] (and developed further by others, including Zerili [10,11] and Moncrief [12]). In the RW gauge, certain 
projections of the MP onto a tensor- harmonic basis, in a specific coordinate system, are taken to vanish. 1 Another 
such "algebraic" gauge proven useful, is the one referred to as the radiation gauge [14], where one sets to zero the 
projection of the MP along a principle null direction of the background black hole geometry. Perturbations of the 
Kerr geometry have been studied almost exclusively through the powerful Tcukolsky formalism [15], in which the 
perturbation is formulated in terms of the Weyl scalars, rather than the metric. A reconstruction procedure for the 
MP, out of the perturbation in the Weyl scalars, has been prescribed by Chrzanowski [14] (with later supplements 
by Wald [16] and Ori [17]). This reconstruction is formulated within the radiation gauge, and relies crucially on this 
choice of gauge. 

Hence, we are in a situation where the singular, "removable" piece of the MP near the particle is given in the 
Lorenz gauge (which best reflects the symmetry of the point singularity), while the full field of the particle can 
only be calculated in the RW or radiation gauges (which bear on the symmetries of the black hole background). 
Much effort has been invested recently in trying to resolve this problematic issue. One approach has been to try 
reformulating the SF, and the subtraction scheme, entirely within the RW gauge (in the Schwarzschild background) 
[18,19]. Another approach incorporated an "intermediate" gauge, in which the MP admits the isotropic Lorenz-gauge 
form near the particle, but gradually approaches the form of the RW (or radiation) gauge as one moves away from 
the particle [8] . Neither approaches have gone very far, partly because of the difficulty in tackling the singular gauge 
transformations involved. The complicated singular nature of these gauge transformations reflects the fact that the 
RW/radiation gauges are far from being natural gauges to describe point particle singularities. For example, is has 
been demonstrated [8] that the radiation-gauge MP develops a singularity along an infinite ray emanating from the 
particle. Indeed, the only actual computation of the gravitational SF carried out so far [20,21] restricted to a special 
case where the RW-gauge MP happens to coincide with the Lorenz-gauge MP (the case of a strictly radial infall into 
a Schwarzschild black hole). 

Of course, the above gauge problem would have never occurred, had we the right tools for calculating the full pertur- 
bation field of the particle in the Lorenz gauge. This would have allowed a direct approach to SF calculations, based 
entirely on the Lorenz gauge and avoiding the above gauge-related complications. In this work we begin developing 
such practical tools for analysis of Lorenz-gauge perturbations. This paper deals with MP of the Schwarzschild black 
hole. It formulates the Lorenz-gauge MP equations in a form accessible to numerical time-evolution treatment, and 
presents calculations of the Lorenz-gauge field in the example of a mass particle moving in a circular geodesic in the 
strong field of the black hole. The paper also discusses, in brief, how we might go about carrying out Lorenz-gauge 
calculations in the Kerr spacetime. This important task is currently subject to intensive study. 

It may be worthwhile to summarize here, in bullet points, the main strengths of our Lorenz-gauge approach to 
black hole perturbations. The following points also give cues to the variety of problems where such an approach can 
prove useful. 

• A regularization scheme for the SF in Kerr spacetime has only been prescribed in the Lorenz gauge. It is not 
clear yet how one goes about formulating and calculating the gravitational SF in other gauges. Having at hand 



1 There is an equivalent formulation of the RW gauge, which does not rely on a tensor-harmonic decomposition [8,13]. In this 
formulation one sets to zero certain combinations of the MP components and their derivatives, in a certain coordinate system. 
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the Lorenz-gauge MP allows one, immediately, to compute the gravitational SF, without having to resort to 
cumbersome gauge adjustments. 

• In our approach one solves directly for the MP components, without resorting to complicated reconstruction 
procedures. This feature becomes important in problems where knowledge of the Weyl scalars (or RW/Zerilli's 
variables, or Moncrief 's variables) is insufficient, and direct access to the MP itself is necessary. 

• Related to the last point is the fact that, in our approach, the MP reconstruction is algebraic (see Sec. HE 
below), and does not involve differentiation of the field variables. This comes to be a great advantage in 
numerical applications, where differentiation often results in loss of numerical accuracy. In particular, in SF 
calculations it is necessary to resolve the MP near the particle with a great precision. The higher the order of 
the derivatives involved in constructing the SF are, the tougher become the resolution requirements. If the SF 
is to be constructed from the Weyl scalars (or from Moncrief 's variables, as in [21]), this construction would 
necessitate taking three successive derivatives — two to reconstruct the MP, and a further one to obtain the force 
exerted by the MP. This has proven to be very demanding computationally [21]. In our approach, one need only 
take a single derivative of the numerical field variables. 

• Yet another advantage of working with the Lorenz-gauge MP components as field variables, is that these behave 
more regularly near point particles than do Tcukolsky's or Moncrief 's variables. This has a simple manifestation 
when considering the multipolc decomposition of the MP: The individual multipole modes of the Lorenz-gauge 
MP arc continuous at the particle; only their first derivatives are discontinues there. The multipole modes of 
Tcukolsky's or Moncrief's variables, on the other hand, are themselves discontinuous at the particle, and so are, 
in general, the modes of the MP in the RW gauge. Obviously, this better regularity of the Lorenz-gauge MP 
comes to be a great advantage when it comes to numerical implementation. 

• The Lorenz-gauge perturbation equations take a fully hyperbolic form. (Compare with the situation in the RW 
or radiation gauges, where the set of perturbation equations split into a subset of hyperbolic field equations, and 
a subset of elliptic equations that constrain the evolution.) This makes the Lorenz-gauge formulation especially 
convenient for numerical applications which are based on time-domain evolution. The supplementary gauge 
conditions indeed take the form of elliptic "constraint" equations, but these can be made to hold automatically. 
We shall discuss this point in Sec. II. 

• The Lorenz gauge condition does not interfere with the local isotropic symmetry of point particle singularities. 
This makes it well suited — in gravitational perturbation theory as in electromagnetism — for study of the local 
field near such particles. Other gauge choices (the RW and radiation gauges are examples) may artificially 
distort this symmetry, requiring a more cautious treatment. 

• The use of (generalized) "harmonic coordinates" , closely related to the choice of Lorenz gauge in a perturbative 
context, has lately been a popular trend in Numerical Relativity (e.g., [22]). It has been partially responsible for 
the significant progress made recently in fully non- linear simulations of binary black hole mergers [23]. There 
is still much to understand about the underlying mathematics that makes this formulation so successful. The 
simpler realm of perturbation theory can provide a good test bed for these ideas, and help gaining insight into 
some of the mathematical features of hyperbolic formulations. 

The obvious down side of the Lorenz-gauge formulation is that one can no longer benefit from the kind of algebraic 
simplicity that has made the RW/radiation gauges so popular. There is no known way to fully decouple between the 
various tensorial components of the Lorenz-gauge MP, and one has to treat the perturbation equations as a coupled 
set. This, seemingly, has discouraged the development of a Lorenz-gauge MP formulation in the past. Part of our 
aim here is to demonstrate that such formulation is tractable, and accessible to numerical implementation, despite 
the lack of algebraic simplicity. 

This paper is structured as follows. In Sec. II we formulate MP theory in the Lorenz gauge. Specializing to the 
case of a point mass source in Schwarzschild geometry, we decompose the field equations (and the gauge conditions) 
into tensor harmonic multipolcs, and obtain a set of 10 hyperbolic wave equations for each multipole mode (each 
l 7 m). This set couples between the various tensorial components of the MP; however, the MP variables are chosen 
in such a way that the equations are decoupled at their principle parts. Using the supplementary gauge conditions 
we modify the original wave equations to incorporate a "constraint damping" scheme: A mechanism designed to 
guarantee that violations of the gauge conditions (due, e.g., to numerical errors) are damped automatically during 
numerical evolution. In Sec. Ill we implement the above formulation for the case of a source particle in a circular 
geodesic orbit around the black hole. We give analytic solutions for the monopole mode of the MP, and for all of 
its axially-symmetric, odd-parity modes. To solve for the rest of the modes, we present a numerical code, based on 
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characteristic time-evolution in 1+1 dimensions. Section IV presents a series of validation tests for our code. In 
particular, we use our Lorenz-gauge solutions to extract the fluxes of energy (in each of the multipole modes) radiated 
to infinity in gravitational waves, and compare them with the values given in the literature. Finally, in Sec. V, we 
discuss the applicability of our approach to more general orbits in Schwarzschild, and to orbits around Kerr black 
holes. 



II. FORMULATION 

A. Linearized Einstein equations in the Lorenz gauge 

Let g^v be the metric in a given "background" spacetime, which we assume to be Ricci-flat [i.e., Rafiig^v) = 0]. 
Let then represent a small gravitational perturbation away from g^, produced by a given energy- momentum 
distribution T a p. Linearization of the Einstein equations, G a/ 3(g fll/ + h^ u ) = 8-7rT a| g, in the perturbation h^ v about 
the background g^ v yields the general form [11] 

□/ia/3 - g a fPh + h, aP + 2R^ a v ^ V ~ ~ + 9 a pK^ V = -167rT Q(3 , (1) 

where a semicolon denotes covariant differentiation in the background metric g^ u , □ = ; a ;A is the covariant 
D'Alambertian operator, h = g^h^ is the trace of h^ Vl and indices are raised and lowered using the back- 
ground metric g^ v . Here, and throughout this paper, we follow the conventions of Ref. [24]; hence, the metric 
signature is ( — h++), the connection coefficients are T^ v = ^g Xc '{g all . v + g av ^ — g^v.a)-, the Ricmann tensor is 
R a XliV = Tl V fi - „ + T^r^ - r^r^, the Ricci tensor and scalar are R afj = R^p and R = R a a , and the 
Einstein equations are G Q( g = R a p — \g a (iR — 87rT Q/ 3. We shall use standard geometrized units, with c = G = 1. 
It is convenient to re-express the MP equations (1) in terms of the new variables 

h a fj = h a p - -g a ph. (2) 

[Note h(= Jiff) = —h; hence h a p is referred to as the "trace-reversed" MP.] Imposing the Lorenz gauge condition, 

h a p = 0, (3) 

the MP equations reduce to the compact form 

Uh afi (x) + 2fl'V /3 / V = -167rT a/3 . (4) 

This is a "linear, diagonal second order hyperbolic" system, which admits a well posed initial-value formulation on 
a spacelike Cauchy hypersurface (see, e.g., Theorem 10.1.2 of [25]). Furthermore, if the gauge conditions (3) arc 
satisfied on the initial Cauchy surface, then they are guaranteed to hold everywhere [assuming that Eqs. (4) are 
satisfied everywhere, and that T ct( g ;/3 = 0]. 2 

The main motivation for this work comes from problems where the source of perturbation can be regarded as a 
"point particle" , with a definite trajectory defined on the background metric. (See [1] for a discussion of how the 
notions of a "point particle" and a "trajectory" can be made sense of in a general-relativistic context.) In this case, 
the source term in the MP equations takes the form 

/oo 
(-.9r 1/2 - x£(T)}u a u p dr (particle case), (5) 
-oo 

where fi is the mass of the particle, g is the determinant of x£(t) denotes the particle's worldline (parametrized 
by proper time r), and u a = dx^/dr is a tangent 4- velocity defined along the worldline. 



2 The conditions (3) do not fully specify the gauge: There is a residual gauge freedom within the family of Lorenz gauges, 
h a /3 — ► h al 3 + £ a;( 3 + £/3 ;a , with any £ M satisfying D£ M = 0. It is easy to verify that the form of both (3) and (4) is invariant 
under such gauge transformations. 
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B. Multipole decomposition 



In the rest of this work we restrict our discussion to perturbations of the Schwarzschild black hole spacetime. The 
line element in the background geometry is then given by 

ds 2 = -fdt 2 + f- x dr 2 + r 2 (d0 2 + sin 2 9 dip 2 ), (6) 

where / = 1 — 2M/r, M is the mass of the black hole, the event horizon is at f(r) — 0, and t,r,6,ip arc the 
standard Schwarzschild coordinates, which we adopt throughout this paper. We shall proceed by decomposing the 
MP into multipole harmonics. This will be achieved by projecting h a fj onto a basis of 2nd-rank tensor harmonics, 
defined, in the background Schwarzschild geometry, on 2-spheres t,r = const [26]. The spherical symmetry of the 
background geometry will guarantee that the individual multipole harmonics ("Z,m modes") are eigenvectors of the 
wave operator in Eq. (4), and hence evolve independently. However, the ten tensorial components of (each l,m mode 
of) the perturbation will generally remain coupled. 

We will adopt here the basis of tensor harmonics defined in Appendix A, which represents a slightly modified 
version of the one introduced by Zerilli [10]. We denote our tensor harmonics by Y^ lm (9, ip; r), where I and m are 
the multipole and azimuthal numbers, respectively, i — 1, . . . , 10 labels the ten elements of the tensorial basis, 3 and 
a, f3 are tensorial indices. 4 The harmonics Y^ lm depend only on the Schwarzschild coordinates 9 and ip, except 
for simple multiplicative powers of r and f(r) that are included in their definition (see Appendix A) — the former 
in the interest of dimensional balance, and the latter for regulating the behavior at the event horizon. With this 
definition, the noncoordinatc-basis components Y^) lm = {g aa g l3t3 ) 1 ^ 2 Y^f j lm (no summation over repeated indices) are 

all dimcnsionless, and become r-independent at the limit M/r — > 0. Also, these harmonics are regular at the event 
horizon [f(r) = 0], in the sense that they attain finite (generally nonzero) values there when expressed in regular, 
"horizon penetrating" coordinates (like the Kruskal coordinates). 
The Y^ l p lm, s constitute an orthonormal set, in the sense that 



dn V a *rf v \YW m ]*Y$ l ' m ' = <M«"W (7) 



(for any i,j = l,..., 10), where rj atJ - = diag (l, /, r~ 2 ,r~ 2 sin~ 2 9), an asterisk denotes complex conjugation, and the 
integration is carried over a 2-sphere of constant r and t. The seven harmonics i — 1, ... ,7 constitute a basis for 
all covariant 2nd-rank symmetric tensors of even parity, while the remaining three harmonics i = 8, 9, 10 span all 
such tensors which are of odd parity. Any covariant 2nd-rank symmetric tensor t a p can hence be expanded as t a p = 

E| m E!!i t(i)lra ( I '.<) y i* m . where the time-radial coefficients are given by t^ lm (r,t) = J dnt a f 3 r) a »r)P"[Y$ lm ]* . 
We expand the trace-reversed metric perturbation in the above tensor harmonics as 



^ = ~ r E E * {l)l h {l)lm {r, t) Y$ lm (9, <p; r). (8) 



The coefficients a^ 1 are introduced for the purpose of simplifying the form of Eqs. (18) below, and are given by 

f 1, 2 = 1,2,3,6, 

0,™ = -= x <^ [1(1 + 1)]- 1 / 2 , 2 = 4,5,8,9, (9) 
[ + i = 7,10, 

where 

A = (l + 2)(l- 1). (10) 



3 At I = and / = 1 there are actually fewer than ten independent basis elements: There are only four independent elements 
at I = 0, and eight independent elements at each of the three I = 1 modes — see Sec. II D below for a more detailed discussion. 

4 To avoid confusion we note that, despite the similar notation, the Y^ lm, s adopted here differ from the basis used by one of 
the authors in [27]. 
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The dimensionless scalar fields h^ lm {r,t) will serve as integration variables for the numerical time-evolution of the 
mode-decomposed perturbation equations. For that reason, the above construction is careful in making sure that 
these fields are well behaved (i.e., attain finite, generally nonzero values) both at spacial infinity and along the event 
horizon. The factor 1/r in front of the expansion (8) was inserted to guarantee h^ Lm oc const as r — > oo. The ft,M im 's 

are guaranteed to be regular at the horizon since both the physical perturbation h a p and the harmonics Y^ lm are 
regular there. 

In a similar manner, we expand the energy-momentum tensor as 



'a/3 



10 

J2J2^ i)lm (r,t)Y^ lm (e, V ;r), (11) 



!,m i=l 



with the time-radial coefficients given by 

T« iro (r,t) = J dnr] a ^ v [Y^ lm }*T a/s . (12) 

For the point particle source with energy- momentum given in Eq. (5), this becomes 

T^ lm {r:x v ) = -^u a u^{x p )^ v {x v )Y^ lm *{n v )5{r - r p ) (particle case), (13) 



where, recall, x p (r) denotes the particle's trajectory, and u a its 4- velocity. 

We now wish to obtain equations for the various fields ftW im (r, t) (numbering ten, for each given Z, m), that are fully 
separated with respect to l,m, and are uncoupled with respect to i at their principle part. This is achieved by first 
substituting both expansions (8) and (11) into the field equations (4), and then constructing certain combinations of 
the resulting equations. For example, to obtain an equation for ftW" 1 , one has to add the tt component of Eq. (4) 
to the rr component of that equation multiplied by f 2 . For some of the other h^ lm, s ne has to combine certain 
derivatives of the field equations. The full list of necessary combinations is given in Appendix B. This procedure 
yields a set of 10 equations, where the principle part of the i'th equation involves solely the i'th MP function, ftW m . 
One proceeds by showing that, in each of the 10 equations, the angular dependence on both sides is simply oc Y lrn , 
and then using the orthogonality of the spherical harmonics to separate the equations into multipole modes. 

To write the resulting separated equations in a convenient form, we introduce the standard 2-dimensional scalar-field 
wave operator (including centrifugal potential), 



D 2d = d + - 



(14) 



where /' = df/dr = 2M/r 2 and v,u are the Eddington-Finkelstein null coordinates, defined by v = t + r* and 
u = t — r*, with dric/dr = f . (Throughout this paper it is to be understood that partial derivatives with respect to 
v or u are taken with fixed iioru, respectively, while partial derivatives with respect to t or r are taken with fixed r 
or t, respectively.) The separated field equations then take the form 

n 2 Jh^ lm + M%h^ lm = 47r/i- 1 [r/ ■/o«]T« ,m = S {i)lm (i = 1,. .., 10). (15) 



0) 

*~\ r-\i rr t~\ t~\ n 1 



The quantities M. ?1 are differential operators, of the first order at most, that couple between the various ft^O's 
(with same l,m). The explicit form of the M. (j)'s can be found in Appendix C. As expected, one finds that the 
seven equations for the even-parity modes /jt 1 '---' 7 ) decouple from the remaining 3 equations for the odd-parity modes 
^(8,9,10). We have j^m = o for any i = 1, . . . , 7 with j = 8, 9, 10, and for any i = 8, 9, 10 with j = 1, . . . , 7. Further 
reduction of the system will be achieved in the next step, where we reemploy the gauge conditions. 



C. Gauge conditions and constraint damping 



The field equations (4) are supplemented by the gauge conditions Z a = h a jf = of Eq. (3). In a mode-decomposed 
form, these four conditions read 

H[ m (r,t) = -h { l } - + f (hf + m/r - ftW/r) = 0, (16a) 
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H l 2 m {r,t) = hf - f (hf - hf^j + (1 - AM/r)h^/r- (f/r) (h {1) - fe (5) - 2fh^ = 0, (16b) 



H l 3 m (r, t) = hf - f (hf + 2h^ /r + 1(1 + 1) /r - U 7) /r) = 0, 



(16c) 



H{ m (r, t) = hf - f (hf + 2m It - fc( 10 )/r) = 0. 



(16d) 



[To obtain these equations, insert the expansion (8) into the equations Z t — 0, Z r = 0, (sin 9 Zg),g + (Z v / sm9), v = 0, 
and Zo :ip — Z V: e = 0, respectively; then show that the angular part in all cases is proportional to Y lm and use the 
orthogonality of the spherical harmonics to separate the equations.] These conditions, relating the various ftW's and 
their first derivatives, are to supplement the separated field equations (15). 

One now faces one of the standard problems of Numerical Relativity: the fact that the combined set of 10 field 
equations and 4 gauge conditions (or "constraints" ) is over-determined. Given initial data, the solutions arc determined 
uniquely, in principle, by evolving the 10 field equations; how do we then make sure that the gauge conditions are 
satisfied as well during the evolution? In our case of Lorenz-gauge perturbations, and in the continuum limit, theory 
has it that if the gauge conditions are satisfied on the initial Cauchy surface, then the field equations will preserve 
them throughout the evolution. To see why this is true, simply take the divergence of Eq. (4), which, with the help 
of the contracted Bianchi identities and assuming Ricci flatness of the background geometry and T a ^f — 0, yields 
a simple homogeneous wave equation for the divergence of h a p: OZ a = 0. Thus, in the continuum limit, if h a p 
satisfies a well-posed initial value problem and Z a = on the initial surface, then Z a = is guaranteed throughout 
the evolution. However, in actual numerical time-evolution implementations, gauge condition violations due to finite 
differentiating and/or roundoff errors can grow out of control even provided appropriate initial data. Moreover, in 
most cases it is practically impossible to devise initial data that satisfy the gauge conditions precisely, and at the 
same time are consistent with the field equations. 

Gundlach et al. [28] recently proposed a general scheme for dealing with the above problem (in the wider context 
of fully non-linear Numerical Relativity), which employs the idea of constraint damping. Adapted to our problem of 
Lorenz-gauge perturbations, the idea would be to add to the linearized Einstein equations, Eqs. (4), a term of the 
form —n(t a Zp + tpZ a ), where n is a positive constant and t a is a future-directed timelike vector field. Obviously, 
the new system of 10+4 equations is equivalent (in the continuum limit) to the original system. Also, the added 
term does not alter the principle part of the field equations and hence does not interfere with their neat hyperbolic 
form. However, the evolution equation for Z a now becomes OZ a — n(t a Zp + tpZa)' 13 = 0, which includes a damping 
term. One then expects that, under a range of circumstances [28], Z a would "automatically" damp to zero during 
the evolution (with a timcscale set by k, if t a is taken to be of unit length). 

Here we will be inspired by the above scheme, but will allow ourselves an amount of freedom in executing it: We 
will seek to add terms oc Z a to our field equations, that would assure efficient damping of the constraints, and at 
the same time would lead to simplification in the final form of the decoupled field equations (putting in mind simple 
numerical implementation). The ultimate "justification" for the specific form selected for the added terms would 
come from numerical experiments with circular orbits, described in the following Sections. 

We take t a = —5^, and add to the field equations (4) a term niS^Zp + 8%Z a ), where k = n(r) = f = 2M/r 2 and 

Z a = (Z r , 2Z r , Zg, Z v ) (Schwarzschild components). This choice simplifies the form of the 10 field equations, and 
leads to partial decoupling as we describe below. Our k is not constant, but is expected to suit its purpose as it varies 
slowly, on a scale of the background curvature, which is much larger than the typical scale for constraint violations in 
evolution of perturbations from a point particle. (As our numerical experiments demonstrate, the latter scale tends 
to relate to the radius of curvature associated with the particle — see, for example, Figs. 10 in the next section.) The 
seemingly odd form of Z a (note it has Z r for its t component) has been selected based on numerical experiments with 
circular orbits. It turns out (see Sec. IV B) to yield efficient damping of all 4 constraints Z a = 0, but we shall not 
attempt here to explain that on theoretical ground. 

At the level of the mode-decomposed equations (17), the addition of the above term to the field equations amounts 
to adding (f /2)H 2 at i = 1, 2, adding (f'/4)H 3 at i = 4, 5, and adding (/'/4)_ff4 at i = 8, 9. This brings the separated 
field equations to their final form: 




(17) 
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M%h (f) = \ffhf + 2^(1 - 4M /0 ~ ^ (5) ) - 2J2 t 1 - 6M / r + 12(M/r) 2 ] ^ (3) + ^(6M/r - l)/i (6) , 

(18a) 



M^O-) = i//^(3) + /' (h<n - ~ h y) + L- 2 (h^ m) + l -U'/r) [(1 - 4M/r)m f (hW 2 fh 



M^h® = \ffhf + ^ [1 - 8M/r + 10(M/r) 2 ] /i (3) - ^ - /i (5) - (1 - AM/r)h^ 
M%h® = \f (/#> - - 1 Z(Z + 1) (//r 2 )^ 2 ) - i/'//r [Sh^ + 2m - + 1(1 + l)h^ 



M % hU) = ^ ( x - 4.5M/r)ftW - ^i(Z + 1) - /i (3) ) + ^(1 - 3M/r) (/(/ + l)/i (6) - /i (7) ) 



/U< 6 ) SCi) = L 



m - - (1 - 4M/r) (/-^ (3) + /i (6) ) 



Mf j} h^ = ^ (1 - 4.5M/r) /i< 9 ) - X (i _ 3M/r) ft (10) , 



(6) 

(18b) 
(18c) 

(18d) 
(18e) 
(18f) 
(18g) 
(18h) 
(18i) 
(18j) 



Recall in these equations / = 1 — 2M/r, f = 2M/r 2 , A = (I + 2) (I — 1), <9 r is taken with fixed t, and 9„ is taken with 
fixed u. 



D. Hierarchical structure of the separated held equations 

Following the above manipulations, the five equations for the even parity modes i = 1,3, 5, 6, 7 no longer couple to 
the remaining two equations for i = 2,4. Similarly, in the odd-parity subset, the two equations for i = 9, 10 no longer 
couple to the third equations, for i = 8. One can then solve the set of field equations (17) in a hierarchical manner, 
starting with the 5 even-parity equations for ftC 1 ' 3 ' 5 ' 6 ' 7 ) and the 2 odd-parity equations for hS 9 ' 10 ', and then using 
the solutions as source terms in the equations for h^ 2 ^ (even parity) and (odd parity). This simplification is 
achieved regardless of the form of the source term T a p. Note that the functions h 1 - 2 ' 4 ^ are those constructing the MP 
components h tr , h t g and ht v [sec Eqs. (20) below], which are associated with the shift vector on the surface t = const. 

The multipole sum in Eq. (8) [and in Eq. (11)] contains the two modes I — 0, 1. At these modes there are fewer 
than ten independent tensor-harmonic basis elements. At I = 0, the "vectorial" elements y( 4 > 5 ' 8 > 9 ) and "tensorial" 
elements y( 7 > 10 ) vanish identically, and the monopole MP is then composed, in general, of only the four "scalar" 
elements y(b 2 > 3 > 6 ) (which are all even-parity). The system of field equations (17) thus reduces, at I = 0, to a 
hierarchical set of 3 coupled equations for h^ 1,3 ' 6 \ plus a single equation for h^ 2 \ As to the dipole, I = 1 mode: Here, 
the two "tensorial" elements y( 7 > 10 ) vanish identically, and one is left with a hierarchical set of 4 + 2 equations for 
the even parity modes, and a second hierarchical set of 1 + 1 equations for the odd parity modes. 

Table I summarizes the hierarchical structure of our separated field equations, for the different values of I. 



8 





Even parity 


Odd parity 


1 = 


1,3,6 ^2 [m = 


,0] 






I = 1 


1,3,5,6 -» 2,4 [m = 


±1] 


9 


^8 [m = 0] 


I > 2 


1,3,5,6,7 -» 2,4 [Z + m 


even] 


9,10 


-» 8 [Z + m odd] 



TABLE I. The hierarchical structure of the decoupled field equations (17) [numbers in this table refer to "j" values of the tenso- 
rial-harmonic modes ftW]: The full set of 10 equations first decouple into two subsets of 7 equations (even parity modes, i = 1, . . . , 7), 
and 3 equations (odd parity modes, i = 8, 9, 10). Both even and odd sectors then further reduce into smaller subsets of equations (in a 
hierarchical sense): For modes with I > 2, the even-parity sector reduces to a subset of 5 equations, for i = 1,3,5,6,7, whose solutions 
are then used as source terms in solving for i = 2,4. Similarly, the odd-parity sector decouples into two (hierarchical) subsets of 2 and 1 
equations. In the monopole and dipole cases the system is simpler: At I = the MP is purely even-parity, and one solves two (hierarchical) 
subsets with 3 and 1 equations. At I = 1 one deals with 4 + 2 equations in the even-parity sector, and 1 + 1 equations in the odd-parity 
sector. When the source of the perturbation is an orbiting point particle, we can choose to work in a Schwarzschild coordinate system in 
which the orbit is confined to the equatorial plane (Q = tt/2). In this case, modes with even values of I + m will be of pure even parity, 
while modes with odd values of I + m will be of pure odd parity. The relevant m modes that contribute to each of the entries of the table, 
in the case of a source particle in an equatorial orbit, are indicated in square brackets. 



E. Reconstruction of the metric perturbation 



Finally, it is useful to have at hand explicit formulas for reconstructing the various components of the original metric 



perturbation h a p, given the functions hS % 
we find 



l (r,t). Using h a p — h a p — ^g a ph, together with Eq. (8) and Appendix A, 

(19) 



73' 



1=0 m=—l 



with 



h l ™ = (h^+fh^ Y lm , 



. Ira _ 



Im 

TV 



Y 



rim 



,lm _ 
'hip — 



h 



rO 



h lm 

h 
h 



Tip 

Ira 



Ira 

ipip 



y-l ft (2)yZm 

r 2 (m-fm) 

r/" 1 (hWYft + hWYft) , 
irf' 1 wi6 (h^Yfe - h^Y{ 
r 2 (f- l h m Y lm + h^Yir + h^Yj^) , 
ir 2 sm9 (h^Y^f -h^Y^ 1 ) , 

^ f -l- h (3) Y lm _ - h {7)yhn _ ^(10)yim 



2 • 2 

r sin 



(20) 



where we have omitted the indices Im from h^ lm for brevity, and where Yy£, Yyg, Kp" 1 , and YJ^ are angular 
functions constructed from the spherical harmonics through 



\rlm 



- r V2 



\rlm 

I T1 



1 



1(1 + 1) 
1 



ylm 

,0 ' 



sin 



-10 ylm 



Y, 



T2 



1(1 + 1) 

j^—^ [sin* (sin" 1 01$") e -sm- 2 6Y; 



(21) 
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Note Eq. (19) gives the MP itself, not its trace-reversed counterpart. In applying Eqs. (20) to I = 0, 1, recall 
^vT,V2 — ^tTt2 = an d ^tTt2 = 0- Finally, note that the trace of the MP is simply given by 

h = g afi h af3 = ( M /r) ^(r 1 ^ - h (6) )Y lm . (22) 

lm 

This implies that the function must vanish at the event horizon (where / = 0), since the trace h must be regular 
at the horizon, and the functions and are both finite there by construction. 

The various fields h l £p are complex quantities. However, the sum ^2 m is guaranteed to be real. From the 
definition of the harmonics yW Jm in App. A, and the decoupled field equations (17), one readily verifies the symmetry 
relation 

fh^Y^ lm ) = fj l (i) Y ^ lm )* (23) 

(namely, the product h^Y^ lm is invariant under simultaneous sign-reversal of m and complex conjugation), which 
is valid for each i and any I. Hence, the sum over modes m in the reconstruction formula (19) can be written in the 
form 

E ^-^; =0+2 E Re (^)' ( 24 ) 

m=— I m=l 

which is manifestly real. In practice, this allows a more economic implementation of the reconstruction scheme: One 
need only compute the m > modes to reconstruct the full MP. 

III. IMPLEMENTATION: CASE OF A PARTICLE IN A CIRCULAR GEODESIC ORBIT 

A. Setup 

Consider a pointlike particle of mass /i, in a circular orbit around a Schwarzschild black hole with mass M » fi. 
Neglecting SF effects, the particle traces a geodesic x a = Xp(r), with four velocity u a = dXp/dr. In what follows we 
work in a Schwarzschild coordinate system t, r, 9, ip at which the orbit is confined to the equatorial plane: 

Xp{r) = [*(r), r = const, 9 = tt/2, <^(t)] . (25) 

This geodesic is completely parametrized by the radius r , or, alternatively, by an "angular velocity" 

lu = dip/dt = \Jm/tI. (26) 
The circular geodesic can also be parameterized by the (conserved) specific energy, given by 

£ = -u t = f Q (l-3M/r )- 1/2 , (27) 
where / = 1 — 2M/r . The four velocity of the particle is then given (in Schwarzschild coordinates) by 

U " = (£/./o)[l,0,0,c4 (28) 

Our goal here is to calculate the physical (stationary) MP associated with this orbiting particle, in the Lorenz 
gauge. 

B. Source terms 

From Eqs. (13) and (15) we obtain the source terms for our decoupled filed equations (17). They read 
S (i)lm (r t) _ Ma (,) i(r _ s x f Y] m * (T/2,wt), i = 1-7 (even), 

i [r,t) - iwta d{r r ) x | Y l e m *{n/2,ujt), i - 8—10 (odd), ^ j 



10 



where the coefficients are given by 



f 1 1 


= " (3) = fS/r , 


in\ 


= = = 0, 




= 2ifomu), 




= r u) 2 , 




= r Q uj 2 [l(l + 1) - 2m 2 


a& 


= 2/ow, 




= 2imrou) 2 . 



Note the special case m = 0: For these axially-symmetric modes the field equations for both M 9 -* and M 10 ** are 
sourceless, and (since these two equations do not couple to any of the other /i«'s) one finds 

e - e } = o. (3D 

The entire axially-symmetric odd-parity perturbation is then described by a single function h( 8 \ satisfying a closed- 
form equation. Since, in the circular orbit case, m = modes are static, this equation is in fact an ordinary 
differential equation (ODE), and is readily solvable analytically. Below we construct the analytic solutions for these 
axially-symmetric odd-parity modes. 

In the even-parity sector, both functions and have vanishing source terms at m = 0. Inspecting the field 
equations in their form (15), with Eqs. (Clb) and (Cld), we observe that and couple not only to each other, 
but also to and . However, since this coupling occurs through t derivative terms, and since in our circular-orbit 
case m = modes are static, these coupling terms in fact vanish, and we find 

htlo = j&Lo = 0- (32) 

Hence, the axially-symmetric even-parity part of the perturbation is described by the 5 functions ftC 1 ' 3 ' 5 ' 6 - 7 ), which 
satisfy a coupled set of ODEs. 



C. Analytic solutions for the axially-symmetric, odd-parity modes 



As we explained above, in the circular-orbit case, axially-symmetric (m = 0) odd-parity modes of the MP are 

— /g\ —(8) 

constructed, at each I, from the single function h m=0 , which is t-independent in this case. Denoting h m ' =0 = <j>i(r), 
the field equation for [Eqs. (17) with (18h)] takes the form 



b'l + Vtir)^ = -4/- 2 4 8 L = A x S(r - r ), 



where a prime denotes d/dr, 



Vi(r) = -r 1 (r) 



1(1 + 1) 4M 



(33) 



(34) 



and the coefficient (3i is given by 

l = -32irfo 1 £L>Y% m=o (e = n/2) 

= J 16/ - 1 ^(-l)('- 1 )/ 2 [^(2/ + l)] 1 / 2 /!!/(?-l)!!, Zodd, 
0, I even. 



(35) 



Hence, we have 4>i(r) = for all modes with even values of I, and we need only consider modes with odd I values. 

Two independent homogeneous solutions to Eq. (33) are given, for / > 2 (we will discuss the mode I = 1 separately 
below), by 



i+i 



71=0 



2+1 



(36) 



n=0 
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where 

x = r/(2M) - 1, (37) 

and the coefficients read 

, l(l + l)(l + n-l)\ .'"f k a l n+k 

a - " (i-n + iy.(n + iy.nv K ~ ^ o fcTT (38) 

These solutions have the following asymptotic behavior at the horizon (r — ► 2 M, x — > 0, / — > 0) and at infinity 

(r, a; — ► oo): 



The solution (/> EH is regular (analytic) at the horizon but diverges at r — > oo, whereas the solution (f>f° is regular at 
r —> oo but irregular at the horizon (it vanishes there, but it is non-differentiable). Matching these solutions at the 
particle's location (r = r ), we hence construct a unique regular, continuous solution for the inhomogeneous Eq. (33): 

6,^- 2M\8,x( ^»<TM, r<r a , 

0,> 2 - 2MAA x | 0r(r)0f H (ro)j r > r0j (41) 

where, recall, A = (I + 2){l — 1). In obtaining this solution we have used the fact that the Wronskian W = EH [0°°]' — 
[0 EH ]'0°° must be constant [since the ODE (33) contains no 4>\ term]. Evaluating it at x = one then easily obtains 
W = -b l /(2M) = — (2M A) -1 . 

At I = 1, the function 0;°° of Eq. (36) fails to be a solution of the homogeneous part of Eq. (33) (although 0™ still 
is a solution). Instead, the general homogeneous solution takes the simple form 

0l=i = ar 2 + b/r, (42) 

where a and b are constants. [Note that for I — 1 the effective potential in Eq. (33) reduces to simply V(r) = — 2/r 2 .] 
The coefficients a and b are determined uniquely by requiring regularity at the horizon 5 and at infinity, and imposing 
continuity at r = ro, along with a 'jump' condition for the derivative there: [$ =1 ]r = Pi=i- This yields 

^=i = 4n>fl=ix| J r/ /°f' r S r °' (43) 
^ 3 y \ (ro/r), r > r , v y 

with /3; = i = 16\/37r/ _1 £ w. 

We can now write down explicitly the solution for the axially-symmetric, odd-parity part of the MP h a p itself. From 
the reconstruction equations (19) [with Eqs. (20)], recalling y v ™ = at m = 0, we find that the only non- vanishing 
axially-symmetric, odd-parity components are h tip = h vt , given by 

h tv [m = 0,odd] = - £ £ Binfl^L yjf=° 
odd ; 

= - E oTTTTn ^(r)sin(9Y;^= (e). (44) 

odd* ^ + 



5 There is a subtlety here: The homogeneous solution oc 1/r is regular everywhere for any finite M and ro. However, we 
do wish to require that the solution remains regular even at the limit M — > (taken with fixed r , which is mathematically 
equivalent to the limit ro — ► oo, taken with fixed M). This excludes the solution oc 1/r at r < ro, as it grows unboundedly at 
the horizon in that limit. 
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The lowest multipole contribution to the sum in Eq. (44) is a "conservative" piece coming from the dipole mode, 
1 = 1. It reads 



h tv \m = 0,/ = 1] = -\n<l>i = i sin0r£= 1,m=o 



4 

(r/r ) 2 , r < r , 



= -2/i/ 1 5wro sin 2 6 x j 



c / ^ >> ( 45 ) 
(ro/r), r > r . 



This agrees with the odd-parity dipole solution first obtained by Zerilli [11], which, as pointed out recently by 
Detweiler and Poisson [29], is a Lorenz-gauge solution. [In comparing our solution (45) with Eq. (4.1) of [29], note 
f ~ 1 £u;r l = [Mr /(1 — SM/ro)] 1 / 2 is the particle's specific angular momentum, denoted L in [29].] This piece of the 
MP describes the shift in the angular momentum content of the perturbation across the surface r — r . 

D. Analytic solution for the monopole (I — 0) mode 

The lowest multipole contribution to the MP comes from the monopole, I = m = mode. This "conservative" 
piece of the MP (which is purely even parity) describes the shift in the mass parameter of the perturbation across 
the surface r = ro- As was the case with the m = 0, odd parity modes considered above, at I = 0, too, the field 
equations (17) simplify enough that one can obtain the solution analytically ; with moderate effort. This solution was 
derived recently by Detweiler and Poisson [29] . 6 For the sake of completeness, and since Detweiler and Poisson do 
not actually write down their solution explicitly (as they are interested mainly in the SF exerted by the monopole, 
not the monopole MP itself), we bring this solution here. 

At r < ro, the non- vanishing components of the Lorenz-gauge monopole perturbation are given by 



•-V<-„) = -' 4/M 

=°( r < rn \ = 

r 3f 



h l t J a {r<r,) = —L—P{r), (46a) 



h l -°(r < r ) - -^jQ(r), (46b) 



where 



h l f e °(r < r ) = (sinfl)- 2 /^=°(r < r„) = AfP(r), (46c) 

A=^ T [M-(r -3M)lnf ], (47) 
3Mr / 

P(r) =r 2 + 2Mr + 4M 2 , Q(r) = r 3 - Mr 2 - 2M 2 r + 12M 3 . (48) 
At r > ro, the solutions read 

h l t f(r > r ) = {3r 3 (r - r) + M 2 (r 2 - 12Mr + 8M 2 ) + 

(r - 3M ) [-rM(r + AM) + rP(r)f In / + 8M 3 ln(r /r)] } , (49a) 

h l =°(r > r ) = - 4 2 ^ to {-^ 3 ^o - 2Mr(r 2 - 6Mr - 10M 2 ) + 3M 2 (r 2 - 12Mr + 8M 2 ) + 

(r - 3M) [5Mr 2 + (r/M)Q(r)f In/ - 8Af 2 (2r - 3M) ln(r /r)] } , (49b) 



6 In fact, Detweiler and Poisson did not obtain their solution by directly tackling the Lorenz-gauge equations. Rather, they 
started with the I — solution derived by Zerilli [11] in a non-Lorenz gauge, and then, essentially, solved the gauge transfor- 
mation equations that take Zerilli's solution to the Lorenz gauge. 
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h l t(r > r ) = (sin 9)- 2 h l =°(r > r ) = {3r 2 M - 80M 2 r + 156M 3 + 

(r - 3M) [-3r 2 - 12Mr + 3(r/M)P(r)f In / + 44M 2 + 24M 2 ln(r /r)] } . (49c) 

It can be readily verified that these solutions match continuously at r = Tq. It can also be checked that they satisfy both 
the field equations (17) and the gauge conditions (16). For this, note that at I = 0, in the circular orbit case, the only 
non- vanishing ftW's are h^ 1 ' 3 ' 6 ^, and use the relations hj^ = 2y/nfi~ 1 r(h t t+f 2 h rr ), h^ = rf~ 1 (h u — f 2 h rr ), 

and hf} = 4^/tt (/ / 'r)h ee . 

The above monopole solution is regular both at the event horizon and at infinity, taking the asymptotic forms 

h\^ = -lAf + 0{p) ) 

f 2 h l =° = \Af + 0(/ 2 ) } as r - 2M (/ - 0), (50) 
h l e = e ° = UM 2 Af + 0(f) J 



/4 =0 = -^/o _1 (l-Wr) I 

h l =° = ^ }+0(l/r 2 ), asr^oo. (51) 

»- 2 ^l = ^ffo^ 1 - 3M/r ) J 

Detweiler and Poisson show [29] that this is a unique Lorenz-gauge solution which is regular both at the event horizon 
and at infinity: Any gauge transformation of this solution within the class of Lorenz gauges would lead to irregular 
behavior at one (or both) of these asymptotic domains. 7 

Physically, the monopole perturbation of Eqs. (46) and (49) describes a shift in the Schwarzschild mass across 
r = tq. This is most clearly evident from Zcrilli's form of the monopole solution [10] (which differs from the above 
Lorenz-gauge solution only by a gauge transformation [29]), where it is easily seen that the geometry described by 
g a + h l =p is that of a Schwarzschild black hole with mass M at r < r , and that of yet another Schwarzschild black 
hole, with mass M + /i£, at r > r$. 

The above Lorenz-gauge monopole solution is plotted in Fig. 1, for a sample of tq values. 

E. Numerical solutions for the rest of the modes 

For the modes considered so far (those with m = and I = 0, 1, 3, 5, . . .) the field equations simplify enough to 
allow a fully analytic treatment (as least within the simplicity of the circular orbit case). To solve for the rest of the 
modes we will resort to numerical methods, employing the full machinery of the formalism developed in Sec. II. In 
the rest of this section we briefly describe our numerical method, and plot a sample of numerical solutions. In the 
next section we shall present various validation tests for our numerical code. 

We need to integrate, numerically, the set of coupled field equations (17), with the source terms given in Eq. 
(29). Recall that for each given l,m, we are facing two sets of 7 (even parity) and 3 (odd parity) equations, which 
couple in the manner described in Table I (but recall that no coupling occurs at the principal parts of the equations). 
We choose to integrate these equations in the time domain, i.e., without introducing a Fourier decomposition. A 
frequency-domain analysis (of the sort employed many times in the past — see, e.g., [30]) is likely to be more efficient, 
numerically, for studying circular orbits. Our strong motivation in developing a time-domain evolution code stems 
from the fact that such a code is readily extensible to any type of orbit, with arbitrary eccentricity. A time-domain 
code can handle radial plunge trajectories as efficiently as it handles circular orbits. Frequency-domains codes, on the 
other hand, quickly loose their efficiency with growing eccentricity, as the number of Fourier modes one must sum 
over grow rapidly with eccentricity. Orbits with eccentricities greater than ~ 0.7 are essentially intractable with a 
frequency-domain treatment [31]. Another reason to opt for a time evolution approach, especially having in mind SF 
applications, is the following: A frequency domain analysis of the field equations requires careful implementation of 
boundary conditions, which becomes increasingly difficult with increasing I values. The technical reason for this is 
explained, e.g., by Hughes in [32]. With a time domain evolution one can avoid this difficulty, either by expanding 



7 Note, however, the peculiar feature of the above solution, htt — * const[= — 2/i£(ro/o)~ 1 ] at r — > oo, which means that the 
perturbed metric, expressed in Schwarzschild coordinates, does not tend to the Minkowski metric at infinity. (Recall, however, 
that this peculiarity merely relates to the choice of gauge; the underlying perturbed geometry is, of course, asymptotically flat.) 
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FIG. 1. The Lorenz-gauge monopole solution, Eqs. (46) and (49), plotted here for three values of the orbital radius, ro = 4M, 7M, 12M. 
By construction (Detweiler and Poisson, Ref. [29]) this solution is continuous across the orbit, and well behaved both at the event horizon 
and at infinity. It is a unique Lorenz-gauge monopole solution with these properties. (The quantity E is the specific energy parameter, 
denoted elsewhere in this paper by £.) 



the spacial boundaries of the numerical domain such as to dismiss any boundary effects, or by using characteristic 
evolution (as described below) which avoids boundaries altogether. Hence, a time domain approach should allow 
better accessibility to the higher I modes. 

A suitable numerical method, for time evolution of the field equations with a particle source represented by a 
delta function, was first presented by Lousto and Price [33]. It has since been implemented in a variety of a cases: 
A radially-falling scalar charge [34], a radially-falling mass particle [21,35], a mass particle in a circular orbit [36], 
and, lately, a mass particle in eccentric and parabolic orbits [37]. In the scalar field case [34] the source term in 
the decoupled field equations [the scalar-field equivalent of our Eqs. (17)] is simply proportional to a delta function 
5[r — T'p(r)]. However, in the gravitational field case, all above works [21,35-37] employed the Regge-Wheeler-Zerilli- 
Moncrief's formulation, where the source involves also terms proportional to dS[r — r p (T)]/dr. As a consequence, 
solutions to the Regge-Wheeler-Zerilli-Moncrief equations turn out discontinuous across the orbit, which, of course, 
complicates considerably the implementation of the above numerical method. One great advantage of working with 
the Lorenz-gauge MP, is that the source term in the field equations involves only a delta function, no derivatives 
thereof — as in the simple scalar case. The solutions are then continuous everywhere, and the implementation of the 
numerical integration scheme of Ref. [33] becomes must simpler. 

Our numerical code is based on characteristic evolution and uses double-null coordinates v = t + r* and u = t — r* 
(like in [34], but unlike, e.g., in [37]). The numerical domain is a two dimensional fixed-step grid, as illustrated in Fig. 
2. The evolution starts with characteristic data on v — vq and u — uq, where we take vq — r„(ro) and uq = — r*(rn). 
That is, we take t = at the initial instance of the evolution [represented by the vortex (vo^uq)], and take the two 
initial null surfaces to be the ingoing and outgoing light rays emanating from the particle at that instance. The 
location of the event horizon on this grid is approximated by a large value of u (with constant v), while large value 
of v (with constant u) approximates the location of null infinity. The trajectory of the particle is represented by the 
vertical line v — u — 2r»(ro), connecting the lower- and upper-most vertices of the grid. Since our numerical evolution 
variables are continuous at the trajectory, there are no complication involved in taking the trajectory to cut through 
grid points, as we do here; in fact, we found this setup most convenient. 

The numerical evolution is based on the finite difference scheme developed in [33] — as simplified to the case where 
the source term contains no derivatives of a delta function, and generalized to deal with a set of simultaneous equations. 
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FIG. 2. Numerical grid for characteristic evolution of the field equations. See the text for description. 

Since the scheme has been described in detail elsewhere, we shall not elaborate on it here, but rather refer the reader 
to [33] or [35]. The finite difference algorithm is second order convergent. The value of the solution at a grid point 
with coordinates (vi,Ui) is approximated, within an error quadratic in the step size, based on the values at the three 
grid points (vi-i,Ui), and which have been solved for in previous steps. For terms containing 

derivatives d v we also use information from the grid points (1^,1^—2) and (fi-i, Ui-2) (which is necessary to maintain 
second order convergence). The numerical evolution proceeds along successive lines v = const, from vq to large values 
of v ( "null infinity" ) , where along each such line the integration proceeds from u = uq to large values of u (the "event 
horizon" ) . 

Obviously, we do not know how to prescribe exact initial data for our problem. (In a our stationary scenario, 
knowledge of the exact solution at one particular moment would be equivalent to having at hand the solution at 
all times!) This is a standard problem of Numerical Relativity. However, unlike the situation in fully non-linear 
simulations of (say) equal-mass black hole mergers, here we have the luxury of being able to stably run our evolution 
for as long as it takes for any effect of "corrupted" initial data to dissipate off the numerical domain, and for the 
solution to settle down to its correct, "physical" value. To explain the basic idea, consider first the case of a scalar 
field, where the evolution is not constrained. Suppose that we prescribe some initial data on the initial surface (e.g., in 
our characteristic initial data formulation, determine the magnitude of the scalar field along v = vq and u = u$), and 
then numerically integrate the inhomogeneous scalar field equation, with a source particle. Suppose that the solution 
we thus obtain represents, in the continuum limit, an exact solution of the inhomogeneous field equation. This solution 
would then differ from the true, "physical" solution by a homogeneous solution of the field equation. We may hence 
view the numerical solution as a superposition of the "physical" solution and a spurious homogeneous perturbation. 
However, homogeneous (i.e., vacuum) perturbations of black hole spacetimes always decay at late time (the "no hair" 
theorem). One need only wait long enough until the spurious waves radiate away and the "true" solution is revealed. 
On theoretical ground, we expect the spurious waves to die off in time with a power-law tail, where the power index 
depends on the multipole number of the mode in question. In the circular orbit case, experiment shows that one 
practically has to wait for about one orbital period before the spurious waves clear out (see [34,21,37], and also Figs. 
3-7 below). 

The above picture becomes slightly more involved in our case, where the evolving fields are components of the MP, 
which are subject to certain constraints in the form of gauge conditions. The initial data are now no longer freely 
specifiable, but are required to satisfy the gauge conditions — or otherwise our solution would not be guaranteed to 
satisfy the gauge conditions even at late time. Alternatively, one may incorporate in the numerical evolution procedure 
a machinery for damping out constraint violations, like the one described above. 

Here we follow the second strategy: We use a version of the field equations that has constraint damping built 
into it [i.e., above Eqs. (17)], and free ourselves from the need to devise exact initial data. This option has two 
obvious advantages: Firstly, the required exact initial data would depend on the source in question, and will have 
to be developed case by case. On the other hand, once an effective constraint-damping scheme is established for 
circular orbits, it should be equally effective in other cases (say, eccentric orbits). Secondly, and more importantly, a 
successful constraint-damping scheme should take care of both initial constraint violations, and constraint violations 
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from numerical errors. Specifying constraint-obeying initial data, on the other hand, would not guarantee, by itself, 
that constraint violations due to numerical (e.g., finite differentiation) errors remain small. 
As initial data for our numerical evolution we take, most simply, 

h (l)lm (v = v ) = h (l)lm (u = uo) = 0, (52) 

for any I, m and each of i — 1, . . . , 10. 

F. Sample numerical results 

Figures 3-7 show typical numerical solutions for the various functions h^ Lm . For these plots we have taken the 
particle to move in a strong-field circular geodesic orbit at ro = 7M (corresponding to a period of T or b ~ 116M). 
The fields are evolved for an amount of time equivalent to 5 T or b, which we find more than sufficient to allow the 
transient spurious waves die off. The linear resolution for these runs is set at a few grid points per M in both v and 
u, which translates to a few hundreds points per wave cycle at m — 1. With this resolution, the magnitude of the 
various h^ lm, s is resolved at a fractional accuracy better than ~ 10 -4 even near the particle's location, where the 
fields' gradients are largest (cf. Fig. 8 in the next section). With these specifications, a single I, m mode takes of order 
a second to run on a standard laptop. 

The graphs in Figs. 3-6 show the various 

h (i)lm, s for l = 2 and m = 12 Plotted 

are the absolute val- 
ues, \h^ lm \ = {[Re(h^ lm )} 2 + [lm(h^ lm )] 2 } 1/2 . The different Figures present different slice-cuts through the 2- 
dimensional numerical grid: r = const(= 7M), t = const, u = const (i.e., an outgoing ray approaching null infinity), 
and v = const (an incoming ray approaching the event horizon). Note that in Fig. 3 we show the fields h^ lm evaluated 
at the particle's location; recall these fields are continuous and admit definite values there. The last plot, in Fig. 7, 
shows the numerical solutions for I — 10 and m = 9, 10. High-Z calculations are more demanding computationally, 
as the spacial scale of variation of the fields decreases with increasing I. In its present from, our code can handle 
multipolcs up to I <~ 20; this should be sufficient for high precision calculations of the SF, through the mode-sum 
scheme. 

Here are some features to notice when examining Figs. 3-7: (i) The early stage of the fields' evolution is dominated 
by spurious waves associated with the imperfection of the initial data. These arc transients, and die off almost entirely 
within one orbital period of evolution time (the rate of decay of the spurious waves seems roughly independent of 
l,m). Typically, one can safely read off the values of the fields after ~ 2 T or b of evolution time. In performing precision 
calculations (e.g., of the SF) for periodic orbits, it is easy (and advisable) to monitor the level of "contamination" from 
transient waves, by comparing the values of the fields at, say, t = T or b, 2 T or b, 3 T or b, . . .. (ii) As emphasized above, 
the fields ftW' m are all continuous through the particle's location. This is manifested in Figs. 4, 5, 6, and 7. Those 
fields /jW™ whose evolution equations have non-vanishing source terms (in our circular-orbit case, i — 1,3,4,6,7,8) 
have discontinuous spacial derivatives across the particle, as expected on theoretical grounds. Those functions that 
are sourceless (i — 2,5,9, 10) have continuous derivatives there. The values of the fields and their gradients at the 
particle provide sufficient information for calculating the gravitational SF. Our code resolves these values with great 
accuracy, (iii) The various /tM im ' s approach finite values at null infinity and toward the event horizon, as they are 
designed to do by construction. These finite values are generally non-zero, except that h^ 3 ' lm vanishes toward both 
ends. The reason for the vanishing of h^' lm at the horizon has been explained above [see the discussion surrounding 
Eq. (22)]. The reason for the vanishing of this function at r — > oo will be discussed in Sec. IV C below. 

Finally, we remind that, given the h^ Lm, s, the MP itself is constructed algebraically through formula (19). 

IV. CODE VALIDATION 

In this Section we present a series of tests we have performed to check the validity of our numerical evolution code. 
These include (i) a test of numerical convergence, (ii) confirmation that the numerical solutions satisfy the Lorenz 
gauge conditions, and (iii) comparison of the flux of energy radiated in gravitational waves to infinity, as extracted 
from our solutions, with the values obtained using other methods. The second and third of these tests check both our 
new formulation of the MP equations, and the validity of its numerical implementation. 
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FIG. 3. Numerical solutions for the (dimensionless) Lorenz-gauge MP functions / l ( l ) ; = 2 . m = 1 . 2 j evaluated at the particle's location, for 
a particle in a circular orbit at r = tq = 7M. In the Lorenz gauge (unlike in the Regge- Wheeler gauge, for example) the mode-decomposed 
MP is continuous at the particle, and has a definite value there. The early stage of the time evolution is dominated by transient spurious 
waves associated with the imperfection of the initial data. This part of the evolution (which, of course, is to be discarded in interpreting 
the physics content of the numerical results) lasts around one orbital period of evolution time, after which the inherent physical behavior 
is unveiled. Our code evolves separately the real and imaginary parts of the complex functions h^ lm , which are both needed to construct 



the full MP [through formula (19)]; for compactness, we show here only the moduli |/jW' m 
mode I = m = 2 is of even parity, and is constructed from the seven even-parity functions hS 1 
constructed from the three remaining functions, M 8 > 9 > 10 ). 
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A. Numerical convergence 

Our code is second-order convergent in the fields . This is demonstrated in Figs. 8 and 9 for the fields ft/ 4 - 1 
and h( 8 \ (We have chosen these two functions for our demonstration here since the evolution of these specific fields 
couples to that of all other fields hS 1 '.) 



B. Preservation of the Lorenz gauge conditions 

As discussed above, we do not impose the Lorenz gauge conditions actively in the numerical evolution scheme. It 
is therefore important to verify that (or, rather, monitor how well) our numerical solutions indeed satisfy these four 
conditions, given by Hi = H2 = -H3 = H4 = [see Eqs. (16)]. 

To check this, our code contains a subroutine that constructs the functions H (out of the ftW's and their first 
derivatives) along specified slice cuts through the two-dimensional numerical grid. If the gauge conditions are fully 
damped along the specified slice cut, we should expect all of these four quantities to converge to zero with decreasing 
numerical step size (and, since the £Ts are calculated to second-order accuracy, we expect this convergence to be 
quadratic). On the other hand, if there exists a finite-size constraint violation, one or more of the ff's should converge 
(quadratically, again) to a non-zero value. 

Fig. 10 shows the behavior of the four functions H along an outgoing ray at late retarded time (u = 3 T or b), in 
the example of the modes I = m = 2 and I = 2, m = 1. The functions Hi, H2, and H3 are constructed from the 
even-parity MP modes, and are not trivially zero for I — m = 2. The function H4 is constructed from the odd-parity 
MP modes, and is not trivially zero for I = 2, m = 1. We normalize the H's (which have units of 1/distance) by 
(mw) _1 , the typical length scale on which the MP varies (away from the particle). The plots demonstrate that, at 
least in the part of the evolution later than t ~ T or b, all four gauge conditions are satisfied with great accuracy, 
to within numerical error. Similar results are obtained when examining the late time behavior along incoming rays 
(nearing the horizon) or along slices of fixed r. 

We emphasize that the non-zero values of the functions H in Figs. 10 are not associated with constraint violations 
that have not yet been fully damped — these would have tended to vanish rapidly with increasing (advanced) time. 
Rather, the finite-size values of the H's are merely discretization errors (coming both from the calculation of the MP 
fields hy\ and from the construction of the if's), which tend to zero with increasing numerical resolution. Similarly, 
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r/M rJM 

FIG. 4. A different slice cut through the solutions of Fig. 3, this time showing the behavior of the fields fe(v' =a > m=1 > 2 on a t=const 
slice, 2.5 T or b into the evolution. The wavy feature on the left- and right-hand sides (small and large values of r*, respectively) are 
associated with the spurious initial waves propagating inward (toward the black hole) and outward (to infinity), and are to be discarded. 
While all functions fcW m are continuous across the particle (located at r = 7M, r* ~ 8.83M), those functions whose corresponding 
sources are non-zero have discontinues r derivatives there. The fields /i( 2 > 5 > 9 ) im , whose sources S^ 2 ' 5,Q ^ lm vanish, have continuous 

derivatives across the particle. The insets expand the peak area of the plots, for better clarity. The code resolves the gradients of the 
MP fields at the particle with good accuracy (cf. Figs. 8 in the next section). These gradients are needed for SF calculations. Note the 
significant damping in the MP amplitude at r„ < 0. This, presumably, is the effect of the well known potential barrier surrounding the 
black hole. The high-frequency spurious waves, on the other hand, penetrate this barrier with ease. 

the large numerical values of the TJ's near the particle should not be interpreted as indicative of a higher level of 
constraint violation there, since at the particle, too, all functions H converge to zero with increasing resolution. 
(Obviously, larger numerical errors arise in the H's near the particle, since the field gradients are larger there.) 

C. Energy flux at the wave zone 

Further tests of our code could be performed by comparing with other calculations of the Lorenz-gauge MP. The 
only such calculation we are aware of was carried out by Pfenning and Poisson [38] , for a particle in a weak- field 
orbit. However, the analysis in [38] is mainly concerned with the SF acting on the particle, and does not provide 
expressions for the MP itself. Another option is to compare gauge invariant quantities derived from the MP. One 
such quantity is the flux of energy radiated to infinity in gravitational waves. For circular orbits in Schwarzschild, 
the energy flux (as distributed among the I, m modes) was computed previously by Poisson [30] via frequency-domain 
numerical integration of the perturbation equations in the standard Regge- Wheeler gauge. More recently, Martel [37] 
reproduced these fluxes (as well as the fluxes from eccentric orbits) using a time-domain analysis of the perturbation 
equations, still in the Regge- Wheeler gauge. Here we shall compare the energy fluxes constructed from our Lorenz- 
gauge solutions, with those obtained by Poisson and Martel. 

We shall first need an expression for the energy flux at infinity in terms of the Lorenz-gauge MP. Let E°° denote 
the total energy per unit time crossing (outward) the 2-sphere r = const where r is very large. We assume that r is 
large enough that the above 2-sphere resides in the "wave-zone" , where the radius of curvature 1Z is much larger than 
the longest wavelength A of the MP, and the perturbation takes the form of plane gravitational waves. The overdot 
in Eqo , and in the expressions below, may stand for either dt or d u (both taken with fixed r) , since the two are equal 
at the wave zone limit. Note also that at the wave zone we can replace h a /3,r (fixed t) with —h a /3 t t (fixed r), and 
covariant derivatives with ordinary (partial) derivatives, the latter two differing by an amount of 0(\/lZ). 

The flux of energy in the gravitational waves can be obtained from Isaacson's effective energy-momentum tensor, 
T a p, as explained, e.g., in Appendix B of Ref. [37]. However, we must use caution here: The standard expression for 
T a p, as derived in [39] and used extensively in the literature (e.g., [37]), assumes that the MP is given in a gauge 
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FIG. 5. Another slice cut through the solutions of Figs. 3 and 
u = const(= 3 T or b). Once again, the insets expand the peak areas 
times ("null infinity"). The fluxes of energy and angular momentum carried by gravitational waves to infinity are straightforwardly ex- 
tracted from these values (if fact, only \h^ lm \ and \h,( 10 ) lrn \ are needed for this purpose), as we discuss in Sec. IVC. Note how at large 
v we seem to have \h^ \ ~ 0, \h^ \ ~ \ ~ |/i' 5 ^|, and |/i( 8 )| ~ |h( 9 'j. In Sec. IVC we explain how these empirical asymptotic 

relations are predicted theoretically. 



where it is traceless. Our Lorenz gauge MP is generally not traceless, 8 and we will need to generalize the expression 
for the effective energy-momentum tensor to the case ft ^ 0. This can be done by repeating Isaacson's derivation, 
starting at his Eq. (2.4) and going through the averaging procedure described in his analysis — but this time keeping 
track of all the trace terms. The result is 



1 



T — - ( h a/3 h 
— 327r : 



(53) 



where (• • •) denotes an average over a region of spacetime much larger than A. The first term in this expression (which 
is invariant under h a p — * h a p) is the standard Isaacson tensor. The second term is a necessary correction when the 
MP has a non-zero trace. The form (53) is now truly gauge invariant (unlike Isaacson's original formula) and can be 
used to extract the energy flux from the Lorenz-gauge MP. Given the effective energy-momentum tensor, the energy 
flux at infinity is given by [37] 



(54) 



where the integration is carried over the above mentioned 2-sphere. Note that the integration over dQ automatically 
takes care of averaging 7^„ over a scale ~ 1Z, so one can effectively ignore the averaging procedure in Eq. (53) when 
evaluating the flux. 

We next consider the distribution of E^ among the different multipoles. To this end, we substitute for the MP 
in Eq. (53) using the expansion (8). In each of the two quadratic terms we formally replace one of the MP factors 
by its complex conjugate (which is allowed since the MP is a real function). One can sort the resulting terms inside 
the integral J dfl into 3 groups, according to their angular dependence: One group of terms comes with Y l 171 Y lm * , 
the other with Y l g m 'Y 1 ™* + sin" 2 6 Y^ m 'Y 1 ™*, and the third with (D 2 Y l ' m ')(D 2 Y lm *) + sin" 2 e(D 1 Y l ' m ')(D 1 Y lm *). 



Using the orthogonality relations (A4), and replacing all derivatives d r by — d t , one obtains E x = J2i m E 
individual multipolar contributions given by 



T'lTi 



with the 



8 Generally, the supplementary gauge condition h = cannot be imposed in addition to the Lorenz gauge condition (3), unless 
spacetime is globally vacuum — which is not our case here. That our solutions indeed have h ^ is demonstrated in Figs. 3-7, 
recalling the trace is constructed from the h^' 's through Eq. (22) . 



20 



MP along v=const=3 


T . 

orb 








10 


lh (7) l 


— 


8 




f 


6 




|h< 4) (/ 


4 




7 








I 2 









|h (6 »| 




I 


\ ^> 





lh (3 >l 



retarded time u (in 7 J 

v orb' 



1=2, m=2 




2.9 



3.1 




retarded time u (in T J 

v orb' 



FIG. 6. One more slice cut through the solutions of Figs. 3, 4 and 5, showing the behavior along an incoming ray, v = const(= 3 T or b). 
With increasing retarded time, approaching the "event horizon", all fields hW m | settle at constant values. Although hard to tell from these 
plots, these values are generally nonzero (with the exception of h^ lm , which tends to zero at the horizon — see the text for explanation). 
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Here, the subscript oo under the h^'s indicates that these functions are to be evaluated at the wave zone. 

The above expression for E™ simplifies by virtue of the gauge conditions. Consider, for example, the condition 
Hi = [Eq. (16a)]: Neglecting non-derivative terms [which at the wave zone are 0(X/TZ) times smaller than the 
derivative terms], and replacing d r — > — dt, we find h^ t — —h^ t — h^ t . Similarly, from H2 = 0, H3 = 0, and H± = 

(8) 



- (2) 

we find, respectively, h)J t = —h 



i {3) h (i] 

'00, t) ' t oo,t 



= -h 



(5) 



:(3) 



and /i^ t = -hlcj. 



Altogether, we thus have 



\h®\* = \h£\* 



|^| 2 =0, 



\h£\ 2 = \h&\ 2 , 



(56) 



(In our stationary scenario we have Ift-^l 2 = m 2 uj 2 \h^ | 2 , which means that the asymptotic relations (56) hold 
between the functions themselves, not only between their time derivatives. That our numerical solutions satisfy 
these relations at large r is easily visible in Figs. 4 and 5. This provides yet another validity test for our code.) Hence, 
Eq. (55) reduces to the final form 
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Note that the flux carried by each individual mode is extracted from a single ftW function: in the case of 

even-parity modes (even I + m) , or U 1 ^ in the case of odd-parity modes (odd I + m). 

Table II lists the values of E^ n for I = 1, . . . , 5, as computed from our numerical Lorenz-gauge solutions based on 
Eq. (57). For this table we set the particle to move on a circular geodesic at r$ = 7.9456M, to allow comparison with 
Poisson and Martel, who provide results for this radius. The functions and h^ 10 ^ are evaluated at v = 40 T or b and 
u = 3 T or b, corresponding to r ~ 2100AT and t ~ 2500M. It is possible to estimate the error we make in extracting 
the flux at a finite r by looking at the residual values of, e.g., \h^ \ and \h^ \ — \h^\, which should vanish at the 
limit r — > 00. If we assume that the finite-r error in TS 7 > is the same as the largest among the residues in |^ 3 '|, 
l^' 1 ^! — |ft>' 2 ^|) and \h^\ — \h( 5 '\, we obtain from our code that the fractional error in the flux, for even-parity modes, 
is typically of order 0.1-0.2%. In the same way, assuming that the finite-r error in /i( 10 ) is the same as the residue in 
\h^ \ — \h( 9 '\, we expect a fractional error of order 0.01-0.3% for the various odd-parity modes. 
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FIG. 7. The MP functions for Z = 10 and m = 9, 10, on a t=const slice (compare with Fig. 4). Spacial gradients of the MP get 

larger with growing I, but at Z = 10 the particle is still perfectly resolvable (our code performs well for multipole numbers up to ~ 20). 
The seemingly noisy features on both ends of the plots are, in fact, numerically robust: They are transient spurious waves of the type also 
apparent in Figs. 3-6, but now occurring at much higher frequencies. 

Our code reproduces the energy fluxes of Poisson [30] and Martel [37] with very good accuracy. Our values are no 
more than ~ 1% off those of Martel, and no more than 0.07% off those of Poisson. The overall energy flux (summed 
over the first 5 modes) agrees with that computed by Poisson with a mere ~ 5 • 10 -5 relative difference. 

V. SUMMARY AND DISCUSSION OF FUTURE APPLICATIONS 

We advocate here a new approach to the computation of the metric perturbation from a small particle orbiting a 
black hole. The approach incorporates the following "principles" : (a) The particle is represented as a delta function 
source term (rather than an extended distribution) in the field equations; (b) The perturbation equations are for- 
mulated and integrated for the MP itself (rather than generating functions thereof); (c) The MP is solved for in the 
Lorenz gauge; (d) The numerical integration of the perturbation equations is carried out in the time domain. 

The main advantages offered by this approach are: (i) In applications requiring knowledge of the MP itself, one 
avoids the difficult issue of MP reconstruction; (ii) The Lorenz-gauge MP provides the most natural description of the 
field near the particle, and is the most convenient to work with in analyzing the singular structure of this field; (iii) 
The Lorenz-gauge MP can be incorporated immediately into existing schemes for calculating the gravitational SF. (iv) 
Time-domain numerical integration is more efficient in dealing with eccentric orbits than traditional frequency-domain 
methods, and is more easily generalizable from one set of orbits to another. 

Here we developed the above approach, and explored its feasibility, for orbits in Schwarzschild. The core of our 
formulation includes the decoupled field equations, Eqs. (17), along with the supplementary gauge conditions, Eqs. 
(16). We implemented this formulation numerically for the case of circular orbits, using a characteristic time-evolution 
code that incorporates a constraint damping scheme. We demonstrated that such a code can efficiently resolve the 
singular field near the particle, in a computationally inexpensive manner. 

Since our code integrates the perturbation equations in the time domain, it is possible to use it for analyzing any 
orbit in Schwarzschild: For any given trajectory x p (t) (which, if geodesic, can always be taken to be equatorial) 
one merely has to specify appropriate source functions S l ^^ m [r p (r)] in Eq. (17), and feed these functions into the 
subroutine in our code that tells the numerical integrator where the particle is at each step of the evolution. The 
numerical properties of the evolution code (stability, convergence rate, constraint damping) are insensitive to the 
choice of source. Since our code performs well with circular orbits, it will perform well with any other orbit. 

Much more challenging is the extension to the Kerr case, which, however, also provides the main motivation for 
the above approach. On a Kerr background, the perturbation equations are no longer separable into multipole modes 
in the time domain. A time-domain integrator will therefore have to evolve the equations in 2+1 dimensions (two 
spacial dimensions, and time; the azimuthal direction remains separable even in the time domain). There exist a few 
working codes for integrating Teukolsky's equation in the time domain, for vacuum perturbations in Kerr [see, e.g., 
[40,41])]. A similar code will have to be developed for integrating the Lorenz-gauge field equations (4), with a particle 
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FIG. 8. Illustration of numerical convergence: Shown here are numerical solutions for (the absolute values of) the functions ^( 4 ) i=2 > m =2 
(left panel) and ^(8)J=2,m=l (j-jg^t panel) along an outgoing ray u = const = 3 T^j, — the solutions also shown in Fig. 5. In each of the 
figures we have superposed 4 numerical solutions, obtained with different numerical resolutions: 4, 8, 16, and 32 grid points per M (this 
corresponds to the linear resolution; the numbers of grid points per unit grid area M 2 are the squares of these values). The insets show, 
greatly magnified, two details from each plot: one near the peak (the particle's location), and the other at the far right end ("null infinity"). 
It clearly appears that the solutions are numerically convergent, and that the convergence rate is faster than linear. Note also that, with 
the range of resolutions used here, the magnitude of the fields is already resolvable to within fractional errors of ~ 1CT 5 . This is the case- 
even near the particle, where the fields' gradients are largest. 



source. The challenge here is not in the expected high computational expense insomuch as it is in the difficulty of 
treating the particle singularity: In 2+1 dimensions the Lorenz gauge perturbation field is no longer continuous and 
finite near the particle (as it is, conveniently, in 1+1 dimensions). Rather, it diverges toward the particle. Although 
this divergence is slow (logarithmic in proper distance from the particle), it poses a serious problem when it comes to 
numerical implementation. 

One possible way around this problem is to implement a procedure reminiscent of the "puncture" scheme often in 
use in Numerical Relativity. To sketch the basic idea, let us represent the field equation (4), symbolically, by 'dh = 6\ 
where '<5' represents a given fixed trajectory on the black hole background, and the left-hand side is understood to 
contain also the Riemann term. We introduce a worldtube around the particle, of some fixed radius p that we keep as 
a control parameter, p should be much smaller than the background's radius of curvature, and much larger than the 
radius of curvature associated with the particle; say, p ~ y/Mfi. The numerical time-evolution integrator than utilizes 
the following algorithm: At any point (in the 2+1-dimensional numerical grid space) outside the above worldtube, the 
integrator evolves the homogeneous equation Oh — as it is. However, at points inside the tube, the integrator solves 
for a different, "punctured" field h punct = h — h s - ms , where h S i ng is an analytic approximation for the field h near the 
singularity, chosen such that ft. pU nct is continuous everywhere inside the worldtube, including on the worldlinc itself. 
Thus, inside the worldtube one solves the equation Oh punct = S — Oh s i ns , where the source term on the right-hand side 
is now distributed, but 'more regular' than the original source. The evolution across the boundary of the worldtube 
proceeds by matching the external solution h to h punct + h s i ng . That this (simple) idea could work in practice has 
been demonstrated recently with a toy model of a scalar field in Schwarzschild [42] . 

As explained in the introduction, our main drive in developing the computational approach of this paper is the 
problem of calculating the SF. With the new computational tools at hand, freeing us from gauge complexities and 
issues of MP reconstruction, we can straightforwardly implement the mode-sum scheme [7]. The scheme, we remind, 
requires the values of each multipole of the Lorenz-gauge MP and its derivatives at the particle. The modes are then 
"regularized" individually, by subtracting a certain quantity, given analytically (the "regularization parameters" [7]), 
from each. The sum over regularized modes then gives the physical SF experienced by the particle. Our current code 
can be used to calculate the SF in this manner, for any orbit in Schwarzschild. One minor technical matter still to be 
addressed is the fact that the regularization parameters prescribed in the literature correspond (somewhat awkwardly) 
to a sca/ar-harmonic decomposition of the force's components, whereas the entities we calculate numerically are tensor 
harmonic components of the MP, yielding the rector-harmonic components of the force. It will be required to obtain 
the appropriate regularization parameters for the vector-harmonic components, or, alternatively, expand the vector 
harmonic components in scalar harmonics in order to allow them to communicate with the standard parameters. 
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FIG. 9. Demonstration of second-order numerical convergence. These four plots provide a quantitative test of the solution's numerical 
convergence rate. The left and right panels of the upper (lower) row correspond, respectively, to the left and right insets in the left (right) 
panels of Fig. 8. The red (pale) line in each of the plots shows the ratio ( | ftW [4pts/M] - ftW [8pts/M] | ) / ( | [8pts/Af] - [16pts/M] j) , 
where i = 4, 8 and ftM [4pts/M], for example, represents the numerical solution obtained with resolution of 4 grid points per M. The blue 
(dark) lines similarly show the ratio ( | hW [8pts/M] - hW [16pts/M] |) / ( | hW [16pts/M] - hW [32pts/iW] |) . Both ratios approach the value 
of 4 with growing resolution, which indicates second-order convergence. 



Either of these options is essentially straightforward to implement. 
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FIG. 10. Verifying that our numerical solutions satisfy the Lorenz gauge conditions. Plotted here are the (absolute values of the) four 
functions Hi, H2, H3, and H4 of Eqs. (16), normalized by (mu)" 1 , along an outgoing ray u = 3 T or i-,. As in previous figures, the particle 
is moving on a circular geodesic orbit with ro = 7M. We consider here the modes I = m = 2 (for Hi, H2, and H3, which involve the 
even-parity MP modes) and I = 2, m = 1 (for H4, which involves the odd-parity MP modes). The Lorenz gauge condition is satisfied if 
(and only if) all four functions H vanish at the limit of vanishing step size. In each of the plots we show, superposed, a sequence of four 
lines, obtained with decreasing step sizes. The insets magnify specific areas in each of the plots. All four functions H appear to converge 
to values very close to zero, indicating that the gauge conditions are satisfied with good accuracy. See the text for further discussion. 



APPENDIX A: BASIS OF TENSOR HARMONICS 



The tensor harmonics Y^ lm adopted in this paper are given, in Schwarzschild coordinates t, r, 6, (p, by 
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TABLE II. Comparison of energy fluxes per Z,m-mode at infinity, as extracted from our Lorenz-gauge MP, with the values obtained 
using other methods. The values in the table give [in units of {fi/M) 2 } for tq = 7.9456, with m < modes folded over onto m > 0. 

The results from Poisson's and Martel's computations are quoted from Ref. [37]. Both Poisson and Martcl extracted their fluxes from the 
MP in the Rcgge- Wheeler gauge, the former using frequency domain analysis and the latter integrating the field equations in the time 
domain. The values under "this paper" are derived using Eq. (57), with the Lorenz-gauge functions Tv-" 7 ) (for modes with even I + m) and 
h( 10 ) (for modes with odd I + m) extracted at v = 40 T or b and u = 3 T OT i, (corresponding to r ~ 2100M and t ~ 2500M). The values in 
square brackets under "Poisson" and "Martel" give the relative difference, in percents, between their results and ours. 
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where / = (1 — 2M/r), Y lm (8, ip) are the standard scalar spherical harmonics, s = sin 6, A = (Z — 1)(Z + 2), and the 
angular operators D\ and D 2 are given by 

D 1 = 2{d e -cot6)d Vl D 2 = d B g- cotO d e - s- 2 d w . (A3) 

The radial factors involving r and / are introduced for dimensional balance and for settling the horizon behavior. 
The harmonics Y^ lm constitute an orthonormal set, in the sense expressed in Eq. (7). This can be readily verified 
based on the identities 
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where the integration is carried over a 2-sphere r=const, an asterisk denotes complex conjugation, and 8 nn i is the 
Kronecker delta. 

APPENDIX B: SEPARATION OF THE FIELD EQUATIONS 

The field equations (15), which arc fully separated with respect to l,m, and uncoupled at their principle part with 
respect to i, are obtained by substituting both expansions (8) and (12) into Eq. (4), and then considering certain 
combinations of the resulting equations and derivatives thereof. The necessary combinations are given in table III. 
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TABLE III. Component combinations used in separating the perturbation equations. Here s = sin 8, and the angular differential 
operators Di and D2 are those defined in Eqs. (A3). Curely brackets represent Schwarzschild components of Eq. (4); thus, "{tt} + f 2 {rr}" , 
for example, stands for "take the tt component of Eq. (4), and add to it the rr component of that equation, multiplied by f 2 ." 



APPENDIX C: MODE-DECOMPOSED LINEARIZED EINSTEIN EQUATIONS IN THE LORENZ 

GAUGE 



We give here the explicit expressions for the terms A4^.n appearing in the original form of the separated field 
equations, Eqs. (15). For the numerical evolution in this paper we use a different form of the separated equations 
[i.e., Eqs. (17), with Eqs. (18), which incorporate a constraint damping scheme], and so the original functions 
are not needed in our analysis. We nevertheless give them here, as they might be useful as a starting point for anyone 
wishing to implement a different numerical scheme. The terms are given by 
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